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Abstract 

In  most  books  the  Delaunay  and  Lagrange  equations  for  the  orbital  elements  are  de¬ 
rived  by  the  Hamilton- Jacobi  method:  one  begins  with  the  two-body  Hamilton  equations 
in  spherical  coordinates,  performs  a  canonical  transformation  to  the  orbital  elements,  and 
obtains  the  Delaunay  system.  A  standard  trick  is  then  used  to  generalise  the  approach 
to  the  N-body  case.  We  re-examine  this  step  and  demonstrate  that  it  contains  an  im¬ 
plicit  condition  which  restricts  the  dynamics  to  a  9(N-l)-dimensional  submanifold  of  the 
12(N-l)-dimensional  space  spanned  by  the  elements  and  their  time  derivatives.  The  tacit 
condition  is  equivalent  to  the  constraint  that  Lagrange  imposed  ”by  hand”  to  remove  the 
excessive  freedom,  when  he  was  deriving  his  system  of  equations  by  variation  of  parame¬ 
ters.  It  is  the  condition  of  the  orbital  elements  being  osculating,  i.e.,  of  the  instantaneous 
ellipse  (or  hyperbola)  being  always  tangential  to  the  physical  velocity.  Imposure  of  any 
supplementary  condition  different  from  the  Lagrange  constraint  (but  compatible  with  the 
equations  of  motion)  is  legitimate  and  will  not  alter  the  physical  trajectory  or  velocity 
(though  will  alter  the  mathematical  form  of  the  planetary  equations). 
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This  freedom  of  nomination  of  the  supplementary  constraint  reveals  a  gauge-type 
internal  symmetry  instilled  into  the  equations  of  celestial  mechanics.  Existence  of  this 
internal  symmetry  has  consequences  for  the  stability  of  numerical  integrators.  Another 
important  aspect  of  this  freedom  is  that  any  gauge  different  from  that  of  Lagrange  makes 
the  Delaunay  system  non-canonical.  In  a  more  general  setting,  when  the  disturbance  de¬ 
pends  not  only  upon  positions  but  also  upon  velocities,  there  is  a  ’’generalised  Lagrange 
gauge”  wherein  the  Delaunay  system  is  symplectic.  This  special  gauge  renders  orbital 
elements  that  are  osculating  in  the  phase  space.  It  coincides  with  the  regular  Lagrange 
gauge  when  the  perturbation  is  velocity-independent. 


1  Euler  and  Lagrange 

1.1  The  history 

The  planetary  equations,  which  describe  the  evolution  of  the  orbital  elements,  constitute  the 
cornerstone  of  the  celestial  mechanics.  Description  of  orbits  in  the  language  of  Keplerian 
elements  (rather  than  in  terms  of  the  Cartesian  coordinates)  is  not  only  physically  illustrative 
but  also  provides  the  sole  means  for  analysis  of  resonant  interactions.  These  equations  exist 
in  a  variety  of  equivalent  forms  (those  of  Lagrange,  Delaunay,  Gauss,  Poincare)  and  can  be 
derived  by  several  different  methods. 

The  earliest  sketch  of  the  method  dates  back  to  Euler’s  paper  of  1748,  which  addresses 
the  perturbations  exerted  upon  one  another  by  Saturn  and  Jupiter.  In  the  publication  on  the 
Lunar  motion,  dated  by  1753,  Euler  derived  the  equations  for  the  longitude  of  the  node,  D  , 
the  inclination,  i  ,  and  the  quantity  p  =  a  (1  —  e^)  .  Time  derivatives  of  these  three  elements 
were  expressed  through  the  components  of  the  disturbing  force.  Sixty  years  later  the  method 
was  amended  by  Gauss  who  wrote  down  similar  equations  for  the  other  three  elements  and, 
thus,  obtained  what  we  now  call  the  Gauss  system  of  planetary  equations.  (The  history  of 
this  scientihc  endeavour  was  studied  by  Subbotin  (1958),  who  insists  that  the  Gauss  system 
of  planetary  equations  should  rather  be  called  Euler  system.).  A  modern  but  still  elementary 
derivation  of  this  system  belongs  to  Burns  (1976). 

In  his  memoires  of  1778,  which  received  an  honourable  prize  from  the  Academie  des  Sciences 
of  Paris,  Lagrange  employed  the  method  of  variation  of  parameters  (VOP)  to  express  the  time 
derivatives  of  the  orbital  elements  through  the  disturbing  functions’  partial  derivatives  with 
respect  to  the  Cartesian  coordinates.  In  his  memo  ire  of  1783  Lagrange  furthered  this  approach, 
while  in  Lagrange  (1808,  1809,  1810)  these  equations  acquired  their  hnal,  closed,  shape:  they 
expressed  the  orbital  elements’  evolution  in  terms  of  the  disturbing  potentials’  derivatives 
with  respect  to  the  elements.  Lagrange’s  derivation  rested  upon  an  explicit  imposure  of  the 
osculation  condition,  i.e.,  of  a  supplementary  vector  equation  (called  the  Lagrange  constraint) 
which  guaranteed  that  the  instantaneous  ellipses  (in  the  case  of  bound  motions)  or  hyperbolae 
(in  the  case  of  flyby  ones)  are  always  tangential  to  the  physical  trajectory.  Though  it  has 
been  long  known  (and,  very  possibly,  appreciated  by  Lagrange  himself)  that  the  choice  of  the 
supplementary  conditions  is  essentially  arbitrary,  and  though  a  couple  of  particular  examples 
of  use  of  non-osculating  elements  appeared  in  the  literature  (Goldreich  1965;  Brumberg  et  al 
1971;  Borderies  &  Longaretti  1987),  a  comprehensive  study  of  the  associated  freedom  has  not 
appeared  until  very  recently  (Efroimsky  2002,  2003). 

In  the  middle  of  the  19th  century  Jacobi  applied  a  canonical-transformation-based  procedure 
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(presently  known  as  the  Hamilton- Jacobi  approach)  to  the  orbital  dynamics,  and  offered  a 
method  of  deriving  the  Lagrange  system.  This  technique  is  currently  regarded  standard  and  is 
offered  in  many  books.  Though  the  mathematical  correctness  of  the  Hamilton-Jacobi  method 
is  beyond  doubt,  its  application  to  celestial  mechanics  contains  an  aspect  that  has  long  been 
overlooked  (at  least,  in  the  astronomical  literature).  This  overlooked  question  is  as  follows: 
where  in  the  Hamilton-Jacobi  derivation  of  the  planetary  equations  is  the  Lagrange  constraint 
tacitly  imposed,  and  what  happens  if  we  impose  a  different  constraint?  This  issue  will  be 
addressed  in  our  paper. 


1.2  The  gauge  freedom 

Mathematically,  we  shall  concentrate  on  the  N-body  problem  of  celestial  mechanics,  problem 
that  for  each  body  can  be  set  as 


"  hr  . 
r  -f  4  -  =  AF 


^  2 


(1) 


AF  being  the  disturbing  force  that  vanishes  in  the  (reduced)  two-body  case,  r  being  the 
position  relative  to  the  primary,  and  fi  standing  for  G{mpianet  +  •  A  solution  to  the 

unperturbed  problem  is  a  Keplerian  ellipse  (or  hyperbola) 


r  =  f(C'i,...,C'6,  t) 


(2) 


parametrised  by  six  constants  (which  may  be,  for  example,  the  Kepler  or  Delaunay  elements). 
In  the  framework  of  the  VOP  approach  it  gives  birth  to  the  ansatz 

r  =  f(Ci(t),...,a(t),t)  ,  (3) 


the  ”  constants”  now  being  time-dependent  and  the  functional  form  of  f  remaining  the  same  as 
in  ©•  Substitution  of  ©  into  o  results  in  three  scalar  equations  for  six  independent  functions 
Ci{t)  .  In  order  to  make  the  problem  defined,  Lagrange  applied  three  extra  conditions 


,  dCj  dt 


(4) 


that  are  often  referred  to  as  ’’the  Lagrange  constraint.”  This  constraint  guarantees  osculation, 
i.e.,  that  the  functional  dependence  of  the  perturbed  velocity  upon  the  ’’constants”  is  the  same 
as  that  of  the  unperturbed  one.  This  happens  because  the  physical  velocity  is 


=  i  + 


df  da 


dCi  dt 


(5) 


where  g  stands  for  the  unperturbed  velocity  that  emerged  in  the  two-body  setting.  This 
velocity  is,  by  dehnition,  a  partial  derivative  of  f  with  respect  to  the  last  variable: 


^(c'l, ... ,  a,  t) 


af(Ci,  ...,^6,  t) 
dt 


(6) 


This  choice  of  supplementary  conditions  is  convenient,  but  not  at  all  necessary.  A  choice  of 
any  other  three  scalar  relations  (consistent  with  one  another  and  with  the  equations  of  motion) 
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will  give  the  same  physical  trajectory,  even  though  the  appropriate  solution  for  non-osculating 
variables  Ci  will  differ  from  the  solution  for  osculating  ones. 

Efroimsky  (2002,  2003)  suggested  to  relax  the  Lagrange  condition  and  to  consider 


E 


9f  da 


dCi  dt 


=  $(C'i,...,C'6,t) 


(7) 


$  being  an  arbitrary  function  of  time,  ’’constants”  Ci  and  their  time  derivatives  of  all  orders. 
For  practical  reasons  it  is  convenient  to  restrict  $  to  a  class  of  functions  that  depend  upon  the 
time  and  the  ’’constants”  only.  (The  dependence  upon  derivatives  would  yield  higher-than- 
hrst-order  time  derivatives  of  the  Ci  in  the  subsequent  developments,  which  would  require 
additional  initial  conditions,  beyond  those  on  r  and  r ,  to  be  specihed  in  order  to  close  the 
system.)  Different  choices  of  $  entail  different  forms  of  equations  for  Ci  and,  therefore, 
different  mathematical  solutions  in  terms  of  these  ’’constants.”  A  transition  from  one  such 
solution  to  another  will,  though,  be  a  mere  reparametrisation  of  the  orbit.  The  physical  orbit 
itself  will  remain  invariant.  Such  invariance  of  the  physical  content  of  a  theory  under  its 
mathematical  reparametrisations  is  called  gauge  symmetry.  On  the  one  hand,  it  is  in  a  close 
analogy  with  the  gradient  invariance  of  the  Maxwell  electrodynamics  and  other  held  theories. 
On  the  other  hand,  it  illustrates  some  general  mathematical  structure  emerging  in  the  ODE 
theory  (Newman  &  Efroimsky  2003). 

If  the  Lagrange  gauge  0  is  hxed,  the  parameters  obey  the  equation 


i:  la.  c,] 


dt 


(8) 


[Cn  Cj]  standing  for  the  unperturbed  (i.e.,  dehned  as  in  the  two-body  case)  Lagrange  brackets: 


.  ,  EL  EL  _  EL  EL 

^  “  dCn  dCj  dCj  dCn 


(9) 


To  arrive  at  formula  (0,  one  should,  according  to  Lagrange  (1778,  1783,  1808,  1809),  differenti¬ 
ate  0.  insert  the  outcome  into  0,  and  then  combine  the  result  with  the  Lagrange  constraint 
0.  (In  the  modern  literature,  this  derivation  can  be  found,  for  example,  in  Brouwer  and 
Clemence  (1961),  Efroimsky  (2002,  2003),  Newman  and  Efroimsky  (2003),  Efroimsky  and 
Goldreich  (2003).) 

In  the  simplest  case  the  perturbing  force  depends  only  upon  positions  and  is  conservative: 
AF  =  dR{a)/dr  .  Then  the  right-hand  side  of  0  will  reduce  to  the  partial  derivative 
of  the  disturbing  function  i?(r)  with  respect  to  Cn  ,  whereafter  inversion  of  the  Lagrange- 
bracket  matrix  will  entail  the  Lagrange  system  of  planetary  equations  (for  Q  being  the  Kepler 
elements)  or  the  Delaunay  system  (for  the  parameters  chosen  as  the  Delaunay  elements). 

As  explained  in  Efroimsky  (2003),  in  an  arbitrary  gauge  $  equation  0  will  generalise  to 
its  gauge-invariant  form 


i:  ic„  c,] 


dt 


(9f  .  ^  di  dg  - 

- AF  — - —  - $ 

dCn  dCn  dt  dCn 


(10) 


the  Lagrange  brackets  [Cn  Cj]  being  still  dehned  through  0.  If  we  agree  that  $  is  a 
function  of  both  time  and  the  parameters  Ci  ,  but  not  of  their  derivatives,  then  the  right-hand 
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side  of  m  will  implicitly  contain  the  first  time  derivatives  of  Ci  .  It  will  then  be  reasonable 
to  move  them  into  the  left-hand  side.  Hence,  m  will  be  reshaped  into 


E 


d?  d4 


dCj 

dt 


c?f  .  -  dl  dg  - 

- _ - _  - $ 

dCn  dCn  dt  dCn 


(11) 


This  is  the  general  form  of  the  gauge-invariant  perturbation  equations  of  celestial  mechanics, 
that  follows  from  the  VOP  method,  for  an  arbitrary  disturbing  force  AF(r,  r,  t)  and  under 
the  simplifying  assumption  that  the  arbitrary  gauge  function  $  is  chosen  to  depend  on  the 
time  and  the  parameters  Ci  ,  but  not  on  their  derivatives. 

For  performing  further  algebraic  developments  of  m  and  (ED,  let  us  decide  what  sort  of 
interactions  will  fall  within  the  realm  of  our  interest.  On  general  grounds,  it  would  be  desirable 
to  deal  with  forces  that  permit  description  in  the  language  of  Lagrangians  and  Hamiltonians. 


2  Delaunay 


2.1  Perturbations  of  Lagrangians  and  Hamiltonians 

Contributions  to  the  disturbing  force  AF  generally  consist  of  two  types,  physical  and  inertial. 
Inputs  can  depend  upon  velocity  as  well  as  upon  positions.  As  motivation  for  this  general¬ 
ization  we  consider  two  practical  examples.  One  is  the  problem  of  orbital  motion  around  a 
processing  planet:  the  orbital  elements  are  dehned  in  the  processing  frame,  while  the  velocity- 
dependent  fictitious  forces  play  the  role  of  the  perturbation  (Goldreich  1965,  Brumberg  et  al 
1971,  Efroimsky  &  Goldreich  2003).  Another  example  is  the  relativistic  two-body  problem 
where  the  relativistic  correction  to  the  force  is  a  function  of  both  velocity  and  position,  as 
explained,  for  example,  in  Brumberg  (1992)  and  Klioner  &  Kopeikin  (1994).  (It  turns  out  that 
in  relativistic  dynamics  even  the  two-body  problem  is  disturbed,  the  relativistic  correction  act¬ 
ing  as  disturbance.  This  yields  the  gauge  symmetry  that  will  cause  ambiguity  in  dehning  the 
orbital  elements  of  a  binary.)  Finally,  we  shall  permit  the  disturbances  to  bear  an  explicit  time 
dependence.  Such  a  level  of  generality  will  enable  us  to  employ  our  formalism  in  noninertial 
coordinate  systems. 

Let  the  unperturbed  Lagrangian  be  r  /2  —  U{r).  The  disturbed  motion  will  be  described 
by  the  new,  perturbed,  Lagrangian 


C  =  — - f/(r)  +  A/:(r,  r,  t)  , 


and  the  appropriately  perturbed  canonical  momentum  and  Hamiltonian: 

dAC 


p  =  r 


P 


dr 


n  =  pr-C  =  ^  +  U  +  An  , 


AH  =  -  AC 


1  fdACV 


2  \  cir  / 

The  Euler-Lagrange  equations  written  for  the  perturbed  Lagrangian  ED  are: 

••  dU 


r  = - +  AF 

ar 


(12) 


(13) 

(14) 


(15) 
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where  the  disturbing  force  is  given  by 


d/\C  _  d  fdAC\ 
df  dt  \  ) 


(16) 


We  see  that  in  the  absence  of  velocity  dependence  the  perturbation  of  the  Lagrangian  plays  the 
role  of  disturbing  function.  Generally,  though,  the  disturbing  force  is  not  equal  to  the  gradient 
of  AC ,  but  has  an  extra  term  called  into  being  by  the  velocity  dependence. 

As  we  already  mentioned,  this  setup  is  sufficiently  generic.  For  example,  it  is  convenient 
for  description  of  a  satellite  orbiting  a  wobbling  planet:  the  inertial  forces,  that  emerge  in  the 
planet-related  noninertial  frame,  will  nicely  £t  in  the  above  formalism. 

It  is  worth  emphasising  once  again  that,  in  the  case  of  velocity-dependent  disturbances, 
the  disturbing  force  is  equal  neither  to  the  gradient  of  the  Lagrangian’s  perturbation  nor  to 
the  gradient  of  negative  Hamiltonian’s  perturbation.  This  is  an  important  thing  to  remember 
when  comparing  results  obtained  by  different  techniques.  For  example,  in  Goldreich  (1965) 
the  word  ”  disturbing  function”  was  used  for  the  negative  perturbation  of  the  Hamiltonian.  For 
this  reason,  the  gradient  of  so  dehned  disturbing  function  was  not  equal  to  the  disturbing  force. 
A  comprehensive  comparison  of  the  currently  developed  theory  with  that  offered  in  Goldreich 
(1965)  will  be  presented  in  a  separate  publication  (Efroimsky  Sz  Goldreich  2003),  where  we 
shall  demonstrate  that  the  method  used  there  was  equivalent  to  hxing  a  special  gauge  (one 
described  below  in  subsection  2.3  of  this  paper). 


2.2  Gauge-invariant  planetary  equations 

Insertion  of  the  generic  force  into  m  will  bring  us: 


i:  [c„  c,] 


dCj 

dt 


df  dAC  _  df  d  dAC\ 

dC„  dr  dC„  ^  gi  ) 


(17) 


If  we  recall  that,  for  a  velocity-dependent  disturbance, 

dAC  dAC  df  dAC  dC  dAC  df  dAC  d{g  +  $) 
^  ^  5F  dCn 

then  equality  m  will  look  like  this: 


(18) 


E  [Cn  C,] 


dCj 

dt 


dAC  _  dAC  d4 
dCn  dr  dCn 


dCn  dt  \  dr  ) 


21  ($  +  m) 

SC'„  t  ^  gf  ) 


.(19) 


After  subsequent  addition  and  subtraction  of  {l/2)d{{d{AC)/dr)^)/dCn  in  the  right-hand 

• 

side,  the  gauge  function  $  will  everywhere  appear  in  the  company  of  -|-  d{AC)  / dr  : 


E 


df  d 


fdAC 
\  dr 


dCj 

dt 


d 


AC  +  - 


fdACV 

\  5r  / 


dg  df  d  d  AC  ^  \  ^  \ 

dCn  dCn  dt  dr  dCn  J  \  dr  ) 


(20) 
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the  sum  in  square  brackets  being  equal  to  —  While  m  expressed  the  VOP  method 
in  the  most  generic  form  it  can  have  in  terms  of  disturbing  forces  AF(r,  r,  t),  equation  (I2()|l 
furnishes  the  most  general  form  in  terms  of  the  Lagrangian  perturbation  A/l(r,  r,  t)  (under 
the  simplifying  assumption  that  the  arbitrary  gauge  function  $  is  set  to  depend  only  upon 
the  time  and  the  parameters  Cj  ,  but  not  upon  their  derivatives). 

The  Lagrange  brackets  in  (nni)  are  gauge- invariant;  they  contain  only  functions  f  and  g 
that  were  dehned  in  the  unperturbed,  two-body,  setting.  This  enables  us  to  exploit  the  well- 
known  expressions  for  the  inverse  of  this  matrix.  These  look  most  simple  (and  are  either  zero 
or  unity)  in  the  case  when  one  chooses  as  the  ’’constants”  the  Delaunay  set  of  orbital  variables. 
As  well  known,  this  simplicity  of  the  Lagrange  and  their  inverse,  Poisson,  brackets  of  the 
Delaunay  elements  is  the  proof  of  these  elements’  canonicity  in  the  unperturbed,  two-body, 
problem.  When  only  a  position-dependent  disturbing  function  R{r)  =  A£(r)  is  ’’turned 
on,”  the  Delaunay  elements  still  remain  canonical,  provided  the  Lagrange  gauge  is  imposed. 
This  happens  because,  as  is  well  known  (Brouwer  &  Clemence  1961),  the  equations  of  motion 
together  with  the  Lagrange  constraint  yield,  in  that  case,  the  following  equation: 

E  |C«  Q]  ^  ^  ,  AC  =  AC  (f(Ci,  ,  Cs ,  ()  )  =  (f(C,,  ... ,  Ce ,  () )  .  (21) 

which,  in  its  turn,  results  in  the  standard  Delaunay  system. 

In  our  case,  though,  the  perturbation  depends  also  upon  velocities;  beside  this,  the  gauge 
$  is  set  arbitrary.  Then  our  equation  (HU)  will  entail  the  gauge-invariant  Lagrange-type 
and  Delaunay- type  systems  of  equations  that  are  presented  in  Appendix  1.  Interestingly,  the 
gauge-invariant  Delaunay-type  system  is,  generally,  non-symplectic.  It  regains  the  canonical 
form  only  in  one  special  gauge  considered  below  (a  gauge  which  coincides  with  the  Lagrange 
gauge  when  the  perturbation  bears  no  velocity  dependence).  This  can  be  proven  either  by 
a  direct  substitution  of  that  special  gauge  condition  into  the  gauge-invariant  Delaunay-type 
system  given  in  Appendix  1.  An  easier  way  would  be  to  £x  the  gauge  already  in  PH).  and  this 
is  what  we  shall  do  in  the  next  subsection. 


2.3  The  generalised  Lagrange  gauge: 

gauge  wherein  the  Delaunay-type  system  becomes  canonical 

We  transformed  (ED  into  (EOD  for  two  reasons:  to  single  out  the  negative  perturbation  of  the 
Hamiltonian  and  to  reveal  the  advantages  of  the  gauge 


$ 


dAL 

dr 


(22) 


which  reduces  to  $  =  0  for  velocity-independent  perturbations.  The  hrst  remarkable  peculiar¬ 
ity  of  is  that  in  this  gauge  the  canonical  momentum  p  is  equal  to  g  (as  can  be  seen  from 
®  and  (O): 

1.-1.  dAC  ^  , 

g  =  r  —  $  =  r-|-  — ^  =  p  .  (23) 

dr 

We  see  that  in  this  gauge  not  the  velocity  but  the  momentum  in  the  disturbed  setting  is  the 
same  function  of  time  and  ’’constants”  as  it  used  to  be  in  the  unperturbed,  two-body,  case. 
Stated  differently,  the  instantaneous  ellipses  (hyperbolae)  defined  in  this  gauge  will  osculate 
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the  orbit  in  the  phase  space.  For  this  reason  our  special  gauge  ()22j)  will  be  called  the 
’’generalised  Lagrange  gauge.” 

Another  good  feature  of  TFh  is  that  in  this  gauge  equation  CT  acquires  an  especially 
simple  form: 


i:  ic„  c,] 


dCj 

dt 


d  AH 

dC„ 


(24) 


whose  advantage  lies  not  only  in  its  brevity,  but  also  in  the  invertibility  of  the  matrix  emerg¬ 
ing  in  its  left-hand  side.  As  already  mentioned  above,  the  gauge  invariance  of  dehnition  Q 
enables  us  to  employ  the  standard  (Lagrange-gauge)  expressions  for  [CnCj\~^  and,  thus,  to 
get  the  planetary  equations  by  inverting  matrix  [CnCj]  in  ()24j).  The  resulting  gauge- invariant 
Lagrange-  and  Delaunay-type  systems  are  presented  in  Appendix  1. 

In  the  special  gauge  (0^.  however,  the  situation  is  much  better.  Comparing  m  with  ()24|1. 
we  see  that  in  the  general  case  of  an  arbitrary  R  =  A£(r,  r,  t)  one  arrives  from  to  the 
same  equations  as  from  m,  except  that  now  they  will  contain  —  ATi  instead  oi  R  =  AC. 
These  will  be  the  Delaunay-type  equation  in  the  generalized  Lagrange  gauge: 

^  _  dAH 
dt  d{  —  Mo) 

dG  _  dAH 
dt  d{  —  u}) 

dH  _  dAH 
dt  d{  —  Q) 

where 

L  =  I'l  _ 

We  see  that  in  this  special  gauge  the  Delaunay-type  equations  indeed  become  canonical,  and 
the  role  of  the  effective  new  Hamiltonian  is  played  exactly  by  the  Hamiltonian  perturbation 
which  emerged  earlier  in  dH- 

Thus  we  have  proven  an  interesting  THEOREM :  Even  though  the  gauge-invariant 
Delaunay-type  system  dm-  EH)  is  not  generally  canonical,  it  becomes  canonical 
in  one  special  gauge  (j22jl  which  we  call  the  ’’generalised  Lagrange  gauge.”  This 
theorem  can  be  proved  in  a  purely  Hamiltonian  language,  as  is  done  in  Section  3.3. 


d{-Mo) 

dAH 

dt 

dL 

d{  —  oj) 

dAH 

dt 

dG 

d{-Q) 

dAH 

dt 

dH 

H  =  a 

(25) 

(26) 

(27) 

(^1  —  cos  i  .  (28) 


3  Hamilton  and  Jacobi 

3.1  The  concept 

A  totally  different  approach  to  derivation  of  the  planetary  equations  is  furnished  by  the  tech¬ 
nique  worked  out  in  1834  -  1835  by  Hamilton  and  rehned  several  years  later  by  Jacobi.  In  the 
lecture  course  shaped  by  1842  and  published  as  a  book  in  1866,  Jacobi  formulated  his  famous 
theorem  and  applied  it  to  the  celestial  motions.  Jacobi  chose  orbital  elements  that  were  some 


combinations  of  the  Keplerian  ones.  His  planetary  equations  can  be  easily  transformed  into 
the  Lagrange  system  by  the  differentiation  chain  rule  (Subbotin  1968).  Later  authors  used  this 
method  for  a  direct  derivation  of  the  Lagrange  and  Delaunay  systems,  and  thus  the  Hamilton- 
Jacobi  approach  became  a  part  and  parcel  of  almost  any  course  in  celestial  mechanics.  To  some 
of  these  sources  we  shall  refer  below.  The  full  list  of  pertinent  references  would  be  endless,  so 
it  is  easier  to  single  out  a  couple  of  books  that  break  the  code  by  offering  alternative  proofs: 
these  exceptions  are  Kaula  (1968)  and  Brouwer  and  Clemence  (1961). 

Brouwer  and  Clemence  (1961)  use  the  VOP  method  (like  in  Lagrange  (1808,  1809,  1810)) 
which  makes  the  imposition  of  the  Lagrange  constraint  explicit.  Kaula  (1968)  undertakes,  by 
means  of  the  differentiation  chain  rule,  a  direct  transition  from  the  Hamilton  equations  in  a 
Cartesian  frame  to  those  in  terms  of  orbital  elements.  As  explained  in  Efroimsky  (2002,  2003), 
in  Kaula’s  treatment  the  Lagrange  constraint  was  imposed  tacitly. 

It  is  far  less  easy  to  understand  where  the  implicit  gauge  fixing  is  used  in  the  Hamilton- 
Jacobi  technique.  This  subtlety  of  the  Hamilton-Jacobi  method  is  so  well  camouflaged  that 
through  the  century  and  a  half  of  the  method’s  life  this  detail  has  never  been  brought  to  light 
(at  least,  in  the  astronomical  literature).  The  necessity  of  such  a  constraint  is  evident:  one  has 
to  choose  one  out  of  inhnitely  many  sets  of  orbital  elements  describing  the  physical  trajectory. 
Typically,  one  prefers  the  set  of  orbital  elements  that  osculates  with  the  trajectory,  so  that  the 
physical  orbit  be  always  tangential  to  the  instantaneous  ellipse,  in  the  case  of  bound  orbits,  or 
to  the  instantaneous  hyperbola,  in  the  case  of  flybys.  This  point  is  most  easily  illustrated  by  the 
following  simple  example  depicted  on  Fig.l.  Consider  two  coplanar  ellipses  with  one  common 
focus.  Let  both  ellipses  rotate,  in  the  same  direction  within  their  plane,  about  the  shared  focus; 
and  let  us  assume  that  the  rotation  of  one  ellipse  is  faster  than  that  of  the  other.  Now  imagine 
that  a  planet  is  located  at  one  of  the  points  of  these  ellipses’  intersection,  and  that  the  rotation 
of  the  ellipses  is  such  that  the  planet  is  always  at  the  instantaneous  point  of  their  intersection. 
One  observer  will  say  that  the  planet  is  swiftly  moving  along  the  slower  rotating  ellipse,  while 
another  observer  will  argue  that  the  planet  is  slowly  moving  along  the  fast-rotating  ellipse. 
Both  viewpoints  are  acceptable,  because  one  can  divide,  in  an  inhnite  number  of  ways,  the 
actual  motion  of  the  planet  into  a  motion  along  some  ellipse  and  a  simultaneous  evolution  of 
that  ellipse.  The  Lagrange  constraint  (ED  singles  out,  of  all  the  multitude  of  evolving  ellipses, 
that  unique  ellipse  which  is  always  tangential  to  the  total  (physical)  velocity  of  the  body.  This 
way  of  gauge  hxing  is  natural  but  not  necessary.  Besides,  as  we  already  mentioned,  the  chosen 
gauge  0  will  not  be  preserved  in  the  course  of  numerical  computations.  (Sometimes  osculating 
elements  do  not  render  an  intuitive  picture  of  the  motion.  In  such  situations  other  elements 
are  preferred.  One  such  example  is  a  circular  orbit  about  an  oblate  planet.  The  osculating 
ellipse  precesses  with  the  angular  velocity  of  the  satellite,  and  its  eccentricity  is  proportional 
to  the  oblateness  factor  J2.  Under  these  circumstances  the  so-called  geometric  elements  are 
more  convenient  than  the  osculating  ones  (Borderies  &  Longaretti  1987).) 

We  remind  the  reader  that  the  Hamilton-Jacobi  treatment  is  based  on  the  simple  facts 
that  the  same  motion  can  be  described  by  different  mutually  interconnected  canonical  sets 
(g,  p,  'H{q,p) )  and  {Q,  P,  T-C*{Q,P) )  ,  and  that  fulhlment  of  the  Hamilton  equations  along 
the  trajectory  makes  the  inhnitesimally  small  quantities 

d9  =  pdq  —  Tidt  (29) 


and 


dO  =  PdQ  -  n*dt 


(30) 
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perfect  differentials.  Subtraction  of  the  former  from  the  latter  shows  that  their  difference, 

-  dw  =  de  -  de  =  PdQ  -  pdq  -  {n*  -  n)  dt  (31) 

is  a  perfect  differential,  too.  Here  q,  p,  Q,  and  P  contain  N  components  each.  If  we  start 
with  a  system  described  by  (g,  p,  'H{q,p) )  ,  it  is  worth  looking  for  such  a  re-parametrisation 
(Q,  P,  T-C*{Q,P))  that  the  new  Hamiltonian  H*  is  constant  in  time,  because  in  these 
variables  the  canonical  equations  simplify.  Especially  convenient  is  to  find  a  transformation 
that  nullihes  the  new  Hamiltonian  Td*  ,  for  in  this  case  the  new  canonical  equations  will  render 


Fig.l.  These  two  coplanar  ellipses  share  one  of  their  foci  and  are  assumed  to  rotate 
about  this  common  focus  in  the  same  direction,  always  remaining  within  their  plane. 
Suppose  that  the  rotation  of  one  ellipse  is  faster  than  that  of  the  other,  and  that  a 
planet  is  located  at  one  of  the  points  of  these  ellipses’  intersection,  P,  and  that  the 
rotation  of  the  ellipses  is  such  that  the  planet  is  always  at  the  instantaneous  point 
of  their  intersection.  We  may  say  that  the  planet  is  swiftly  moving  along  the  slower 
rotating  ellipse,  while  it  would  be  equally  legitimate  to  state  that  the  planet  is  slowly 
moving  along  the  fast-rotating  ellipse.  Both  interpretations  are  valid,  because  one 
can  divide,  in  an  infinite  number  of  ways,  the  actual  motion  of  the  planet  into  a 
motion  along  some  ellipse  and  a  simultaneous  evolution  of  that  ellipse.  The  Lagrange 
constraint  (lU  singles  out,  of  all  the  multitude  of  evolving  ellipses,  that  unique  ellipse 
which  is  always  tangential  to  the  total,  physical,  velocity  of  the  planet. 
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the  variables  {Q,  P)  constant.  One  way  of  seeking  such  transformations  is  to  consider  W 
as  a  function  of  only  q  ,  Q  ,  and  t .  Under  this  assertion,  the  above  equation  will  entail; 


dW 


dt 


whence 


P  =  - 


dW 

W 

dW 


dQ 


dW 

dq 


dq  =  P dQ  —  pdq  +  {Td  —  Ti*)  dt 


P  = 


dW 


n 


dW 


dQ  '  "  dq  '  dt 

The  function  W  can  be  then  found  by  solving  the  Jacobi  equation 


=  n* 


(32) 

(33) 


,  f  9W  ' 


+ 


dW 


=  n* 


(34) 


where  Td*  is  a  constant.  It  is  very  convenient  to  make  it  equal  to  zero.  Then  the  above 
equation  can  be  easily  solved  in  the  unperturbed  (reduced)  two-body  setting.  This  solution, 
which  has  long  been  known,  is  presented,  in  a  very  compact  form,  in  Appendix  2.  It  turns 
out  that,  if  the  spherical  coordinates  and  their  conjugate  momenta  are  taken  as  a  starting 
point,  then  the  eventual  canonical  variables  Q,  P  corresponding  to  Td*{Q,  P)  =  0  are  the 
Delaunay  elements: 


Qi 

=  L  = 

^/JTd  , 

Pi  = 

-Mo 

Q2 

=  G  = 

(1  -  e2) 

P2  = 

—  u 

Q3 

=  H  = 

J pa  {1  —  e^)  cos  i  , 

Pz  = 

-  Q 

(35) 


3.2  Where  can  free  cheese  be  found? 


The  transition  from  two-body  to  N-body  celestial  mechanics  is  presented  in  numerous  books. 
However,  none  of  them  explain  how  the  Lagrange  constraint  is  implicitly  involved  in  the  for¬ 
malism. 

Before  we  move  on,  let  us  cast  a  look  back  at  what  has  been  accomplished  in  the  two-body 
case.  We  started  out  with  a  Hamiltonian  problem  {q,  p,  Td)  and  re-formulated  its  equations 
of  motion 


Q  = 


dn 

dp 


P  = 


m 

dq 


(36) 


in  terms  of  another  set  ( Q,  P,  Td* )  ; 


q  =  0(g,  P,  t) 
P  =  i’iQ,  P,  t) 


(37) 


so  that  the  above  equations  are  mathematically  equivalent  to  the  new  system 


Q 


dn* 

'W  ’ 


p 


dn* 

'W 


(38) 


The  simple  nature  of  the  2-body  setting  enabled  us  to  carry  out  this  transition  so  that  our  new 
Hamiltonian  n*  vanishes  and  the  variables  Q  and  P  are,  therefore,  constants.  This  was 
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achieved  by  means  of  a  transformation-generating  function  W{q,  Q,  t)  obeying  the  Jacobi 
equation  After  formula  P|1  for  this  function  was  written  down,  the  explicit  form  of 

dependence  (PTj)  can  be  found  through  the  relations  P  =  —  dW/dQ  .  This  is  given  by  fj92|l . 

To  make  this  machinery  function  in  an  N-body  setting,  let  us  first  consider  a  disturbed 
two-body  case.  The  number  of  degrees  of  freedom  is  still  the  same  (three  coordinates  q  and 
three  conjugate  momenta  p),  but  the  initial  Hamiltonian  is  perturbed: 


Q 


din  +  A7f) 

dp 


p 


djn  +  AH) 

dq 


(39) 


Trying  to  implement  the  Hamilton- Jacobi  method  -  (jHl)),  for  the  new  Hamiltonians  in  -|- 
/Sn) ,  {n*  +  An)  and  for  some  generating  function  l/(g,  Q,  t)  ,  we  shall  arrive  at 


dV  ,  dV 


p  = 


dq 


dq  =  P  dQ  —  p  dq  +  [{n  +  An) 


p  = 


dq 

dV 


n  +  An 


dV 


^  M  A  ^ 


dq 


dt 


in*  +  An)  ]  dt  , 
=  n*  +  An  , 


(40) 

(41) 

(42) 


We  see  that  V  obeys  the  same  equation  as  W  and,  therefore,  may  be  chosen  to  coincide 
with  it.  Hence,  the  new,  perturbed,  solution  (g,  p)  will  be  expressed  through  the  perturbed 
’’constants”  Qit)  and  P(t)  in  the  same  fashion  as  the  old,  undisturbed,  q  and  p  were 
expressed  through  the  old  constants  Q  and  P  : 


q  =  (/>(Q(f),  Pit),  t)  , 
p  =  'ipiQit),  Pit),  t)  , 


(43) 


0  and  0  being  the  same  functions  as  those  in  (EZD-  Benehtting  from  this  form-invariance,  one 
can  now  master  the  N-particle  problem.  To  this  end,  one  should  choose  the  transformation¬ 
generating  function  V  to  be  additive  over  the  particles,  whereafter  the  content  of  subsection 
3.1  shall  be  repeated  verbatim  for  each  of  the  bodies,  separately.  In  the  end  of  this  endeavour 
one  will  arrive  to  —  1  Delaunay  sets  similar  to  (El,  except  that  now  these  sets  will  be 
constituted  by  instantaneous  orbital  elements.  The  extension  of  the  preceding  subsection’s 
content  to  the  N-body  case  seems  to  be  most  straightforward  and  to  involve  no  additional 
assumptions.  To  dispel  this  illusion,  two  things  should  be  emphasised.  One,  self-evident,  fact 
is  that  the  quantities  Q  and  P  are  no  longer  conserved  after  the  disturbance  An  is  added 
to  the  zero  Hamiltonian  n*  .  The  second  circumstance  is  that  a  change  in  a  Hamiltonian 
implies  an  appropriate  alteration  of  the  Lagrangian.  In  the  simple  case  of  An  being  a 
function  of  the  coordinates  and  time  only  (not  of  velocities  or  momenta),  addition  of  An  to 
the  Hamiltonian  implies  addition  of  its  opposite  to  the  Lagrangian.  Since  this  extra  term  bears 
no  dependence  upon  velocities,  the  expressions  for  momenta  through  the  coordinates  and  time 
will  stay  form-invariant.  Hence  (if  the  Lagrangian  is  not  singular),  the  functional  dependence 
of  the  velocities  upon  the  coordinates  and  momenta  will,  also,  preserve  their  functional  form 
viq,  p,  t)  : 


without  perturbation  : 


^  dCiq,q,t) 

dq 


q  =  viq,p,t)  , 


12 


(44) 


with  perturbation  : 


9  mg,  q,  t)  +  A£(g,  t)) 
dq 


dC{q,  q,  t) 
dq 


q  =  v{q,  p,  t)  , 


where  the  new,  perturbed,  dependence  q  =  v[q{Q{t),  P{t),  t),  p{Q{t),  P{t),  t),  t]  has  the 
same  functional  form  as  the  old  one,  q  =  v[q{Q,  P,  t),  p{Q,  P,  t),  t]  .  Together  with  (HHll. 
this  means  that  the  dependence  of  the  new  q  upon  the  new  P{t)  and  Q{t)  will  have  the 
same  functional  form  as  the  dependence  of  the  old  q  upon  the  constants  Q  and  P  : 

^  q{Q{t),  P{t),t)  =  ^  q{Q{t),  P{t),t)  .  (45) 


In  other  words. 


6 


E 


0 


(46) 


where  Di  denotes  the  set  of  perturbed  variables  (Q(t),  Pit))  .  In  the  astronomical  applica¬ 
tions,  Di  may  stand  for  the  Delaunay  set. 

This  is  the  implicit  condition  under  which  the  Hamilton-Jacobi  method  works  (in  the  above 
case  of  velocity-independent  disturbance).  Violation  of  (HHll  would  invalidate  our  cornerstone 
assumption  (EHD-  This  circumstance  imposes  a  severe  restriction  on  the  applicability  of  the 
Hamilton-Jacobi  theory.  In  the  astronomical  context,  this  means  that  the  Delaunay  elements 
o  must  be  osculating.  Indeed,  if  Di  denote  a  set  of  orbital  elements,  then  expression  (HF)|l 
is  equivalent  to  the  Lagrange  constraint  discussed  in  Section  1.  There  the  constraint  was 
imposed  upon  the  Keplerian  elements;  however,  its  equivalence  to  (sni),  which  is  written  in 
terms  of  the  Delaunay  variables,  can  be  easily  proven  by  the  differentiation  chain  rule. 


3.3  The  case  of  momentum-dependent  disturbances 

When  the  perturbation  of  the  Lagrangian  depends  also  upon  velocities  (and,  therefore,  the 
Hamiltonian  perturbation  carries  dependence  upon  the  canonical  momenta),  the  special  gauge 
(I22j)  wherein  the  Delaunay-type  system  preserves  its  canonicity  differs  from  the  Lagrange  gauge. 
This  was  proven  in  subsection  2.3  in  the  Lagrangian  language.  Now  we  shall  study  this  in 
Hamiltonian  terms.  Our  explanation  will  be  sufficiently  general  and  will  surpass  the  celestial- 
mechanics  setting.  For  this  reason  we  shall  use  notations  q,  p  ,  not  r,  p  .  The  development 
will,  as  ever,  begin  with  an  unperturbed  system  described  by  canonical  variables  obeying 


Q  = 


m 

dp 


p 


dn 

dq 


(47) 


This  dynamics  may  be  re-formulated  in  terms  of  the  new  quantities  {Q,  P)  : 


q  =  0(Q,  P,  t) 


p  =  P,  t) 

so  that  the  Hamiltonian  equations  (P7|)  are  equivalent  to 


Q  = 


dn* 

'W 


p  =  - 


dn* 


(48) 


(49) 


13 


For  simplicity,  we  shall  assume  that  Ti.*  is  zero.  Then  the  new  canonical  variables  will  play 
the  role  of  adjustable  constants  upon  which  the  solution  (pH|l  of  dlZD  depends. 

We  now  wish  to  know  under  what  circumstances  a  modihed  canonical  system 


Q 


d{n  +  AH) 
dp 


P 


will  be  satished  by  the  solution 


din  +  A7f) 
dq 


AH  =  AHiq,  p,  t)  (50) 


q  =  0(Q(t),  Pit),  t) 
p  =  p)iQit),  Pit),  t) 


(51) 


of  the  same  functional  form  as 

Q  = 


but  with  time-dependent  parameters  obeying 

dAU  ■  dAH 


dP 


P  = 


dQ 


(52) 


This  situation  is  of  a  more  general  sort  than  that  addressed  in  subsection  3.2,  in  that  the 
perturbation  An  now  depends  also  upon  the  momentum. 

Under  the  assumption  of  being  (at  least,  locally)  invertible,  substitution  of  the  equal¬ 
ities 


Q  = 


dAn 

dP 


and 


P  = 


dAn 


into  the  expression  for  velocity 


dq 


leads  to 


dq 


'  dq  dq  dq  dq  ' 
~  dQ 


dAn  dq 

dAn  dp 

dp  dP 

(53) 

dq  dP 

dAn  dq 

dAn  dp 

(54) 

dq  dQ 

dp  dQ 

dq 

p 

dP 

(55) 

dAn 

dq 

(  dq  dp  dq  dp  ^ 

dAn 

(56) 

[dQ  dP  dP  dQ ) 

dp 

Here  the  coefficient  accompanying  dAn/dq  identically  vanishes,  while  that  accompanying 
dAn/ dp  coincides  with  the  Jacobian  of  the  canonical  transformation  and  is,  therefore,  unity: 


dq  dp  dq  dp 


So  if  we  introduce,  in  the  spirit  of  Q,  notation 


then  (jSEI)  will  lead  to 


dq 

®  "  a 


(dAH\ 
^  ^  ^  [  ap  ) 


(57) 


(58) 


(59) 


q,t 
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Expression  establishes  the  link  between  the  regular  VOP  method  and  the  canonical  treat¬ 
ment.  It  shows  that,  to  preserve  the  symplectic  description,  one  must  always  choose  a  particular 
gauge  $  =  d  AH /dp  .  Needless  to  say,  this  is  exactly  the  generalised  Lagrange  gauge 
discussed  in  subsection  2.3.  A  direct,  though  very  short,  proof  is  as  follows. 

On  the  one  hand,  the  Hamilton  equation  for  the  perturbed  Hamiltonian  m  is: 


d  (H  +  AH)  dAH 

^  ^  ’’  + 

while,  on  the  other  hand,  the  dehnition  of  momentum  entails,  for  the  Lagrangian  m-- 

d  mq,  q,  t)  +  AC{q,  q,t))  .  ^  d  AC 

p  = - -  =  q  +  . 

dq  dq 

By  comparing  the  latter  with  the  former  we  arrive  at: 


(60) 


(61) 


V  9P  /,,,  V  ),,, 


(62) 


which  coincides  with  (1221).  Thus  we  see  that  transformation  (gHl)  being  canonical  is  equivalent  to 
the  partition  of  the  physical  velocity  g  in  a  manner  prescribed  by  dSHD,  where  ^  =  d  AH  /  dp  . 
This  is  equivalent  to  our  Theorem  from  subsection  2.3.  Evidently,  for  disturbances  dependent 
solely  upon  the  coordinates,  we  return  to  the  case  explained  in  subsection  3.2  (equations  ipHl 
-  l4f)j)L  in  that  case,  the  Hamiltonian  formulation  of  the  problem  demanded  imposition  of  the 
Lagrange  constraint 

To  draw  to  a  close,  the  generalised  Lagrange  constraint,  $  =  —  dAC/dq  ,  is  stiffly 
embedded  in  the  Hamilton-Jacobi  technique.  Hence  this  technique  breaks  the  gauge  invariance 
and  is  unht  (at  least,  in  its  straightforward  form)  to  describe  the  gauge  symmetry  of  the 
planetary  equations.  It  is  necessary  to  sacrihce  gauge  freedom  by  imposing  the  generalised 
Lagrange  constraint  to  make  use  of  the  Hamilton-Jacobi  development. 

In  this  special  gauge,  the  perturbed  momentum  coincides  with  the  unperturbed  one  (which 
was  equal  to  g  ).  Indeed,  we  can  rewrite  (EH)  as 


d  {C{q,  q,  t)  AC{q,  q,  t) ) 
dq 


q  -  ^  =  g 


(63) 


which  means  that,  in  the  astronomical  implementation  of  this  theory,  the  Hamilton-Jacobi 
treatment  necessarily  demands  the  orbital  elements  to  osculate  in  the  phase  space.  Naturally, 
this  demand  reduces  to  that  of  regular  osculation  in  the  simple  case  of  velocity- independent 
AC  that  was  explored  in  subsection  3.2. 


4  Conclusions 

We  have  studied,  in  an  arbitrary  gauge,  the  VOP  method  in  celestial  mechanics  in  the  case 
when  the  perturbation  depends  on  both  positions  and  velocities.  Such  situations  emerge  when 
relativistic  corrections  to  the  Newton  law  are  taken  into  account  or  when  the  VOP  method 
is  employed  in  noninertial  frames  of  reference  (a  satellite  orbiting  a  processing  planet  being 
one  such  example).  The  gauge-invariant  (and  generalised  to  the  case  of  velocity-dependent 
disturbances)  Delaunay-type  system  of  equations  is  not  canonical.  We,  though,  have  proven 
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a  theorem  establishing  a  particular  gauge  (which  coincides  with  the  Lagrange  gauge  in  the 
absence  of  velocity  dependence  of  the  perturbation)  that  renders  this  system  canonical.  We 
called  that  gauge  the  ’’generalised  Lagrange  gauge.” 

We  have  explained  where  the  Lagrange  constraint  tacitly  enters  the  Hamilton-Jacobi  deriva¬ 
tion  of  the  Delaunay  equations.  This  constraint  turns  out  to  be  an  inseparable  (though  not 
easily  visible)  part  of  the  method:  in  the  case  of  momentum-independent  disturbances,  the 
N-body  generalisation  of  the  two-body  Hamilton-Jacobi  technique  is  legitimate  only  if  we  use 
orbital  elements  that  are  osculating,  i.e.,  if  we  exploit  only  the  instantaneous  ellipses  (or  hy¬ 
perbolae,  in  the  flyby  case)  that  are  always  tangential  to  the  velocity  vector.  Oddly  enough, 
an  explicit  mention  of  this  circumstance  has  not  appeared  in  the  astronomical  literature  (at 
least,  to  the  best  of  our  knowledge). 

In  the  case  of  momentum-dependent  disturbances,  the  above  restriction  generalises,  in 
that  the  instantaneous  ellipses  (hyperbolae)  must  be  osculating  in  the  phase  space.  This  is 
equivalent  to  the  imposition  of  the  generalised  Lagrange  gauge. 

Comparing  the  good  old  VOP  method  with  that  based  on  the  Jacobi  theorem,  we  have 
to  acknowledge  that  the  elegance  of  the  latter  does  not  outweigh  the  power  of  the  former.  If 
we  decide  to  explore  the  inhnite  multitude  of  gauges  or  to  study  the  numerical-error-invoked 
gauge  drift,  we  shall  not  be  able  to  employ  the  Hamilton-Jacobi  theory  without  additional 
structure.  However,  the  direct  VOP  method  unencumbered  with  the  canonicity  demand  will 
immediately  yield  gauge-invariant  equations  for  the  Delaunay  elements  obeying  an  arbitrary 
gauge  condition 


(64) 


$  being  some  function  of  time  and  elements  Di  .  In  Efroimsky  (2002)  these  equations  were 
written  down  for  the  case  of  velocity-independent  perturbation.  If  the  disturbing  force  depends 
also  upon  velocities,  the  Delaunay-type  equations  will  acquire  even  more  terms  and  will  read 
as  dZU-IZEI).  In  the  simple  case  of  a  velocity-independent  disturbance,  any  supplementary  con¬ 
dition  different  from  that  of  Lagrange  will  drive  the  Delaunay  system  away  from  its  canonical 
form.  If  we  permit  the  disturbing  force  to  depend  also  upon  velocities,  the  Delaunay  equations 
will  retain  their  canonicity  only  in  the  generalised  Lagrange  gauge. 

In  the  language  of  modern  physics  this  may  be  put  in  the  following  wording.  N-body 
celestial  mechanics  is  a  gauge  theory  but  is  not  genuinely  symplectic  insofar  as  the  language 
of  orbital  elements  is  used.  It,  though,  becomes  canonical  in  the  generalised  Lagrange  gauge. 

The  applications  of  this  formalism  to  motions  in  non-inertial  frames  of  reference  will  be 
studied  in  Efroimsky  &  Goldreich  (2003).  Some  other  applications  were  addressed  in  Slabinski 
(2003). 
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Appendix  1: 

Gauge-invariant  equations  of  Lagrange  and  Delaunay  types 


We  present  the  gauge-invariant  Lagrange-type  equations.  They  follow  from  (HH)  if  we  take  into 
account  the  gauge- invariance  of  matrix  [CiCj]  dehned  by  Q-  We  denote  by  ATi  the  pertur¬ 
bation  of  the  Hamiltonian,  connected  through  m  with  that  of  the  Lagrangian.  The  latter, 
in  its  turn,  is  connected  through  (cni)  with  the  disturbing  force  (and  acts  as  the  customary 
disturbing  function  when  the  perturbations  are  devoid  of  velocity  dependence). 


da  2 

'd{-An) 

dAC 

d  / 

dt  na 

dMo 

dr 

dMo  \ 

(* 

dg 

1 

a?  ) 

dMo 

de 

1-e^ 

ra( 

-AH) 

dAC  d 

+ 

dt 

na^  e 

dMo 

dr  da 

$  + 


dAC 

dr 


_  di  d 
dMo  dt 

dAC\  _ 

dr  ) 


$  + 


dAC\ 
d?  ) 


(65) 


dAC\ 
d?  ) 


(1  - 

no?  e 


dg  d?  ^  d AC\ 
dMo  dMo  dt  \  d?  ) 

d  ( -  AH)  _  dAC  ^  -  dAC\ 
duj  dr  doj  y  dr  j 


(66) 


dAC\ 

dC  J 


du! 


^  A  ($  dAC\ 

du  dt  \  d?  ) 


du 

dt 


—  cos  i 

n  (1  —  sin  i 


'd{-An) 

dAC 

d 

di 

dr 

di 

V  dr  J 

dAC\ 

dr  ) 


di 


^  dAC\ 

di  dt  \  (9F  / 


+ 


(67) 


(1  -  e2)V2 

'd^-AH) 

dAC 

d 

U  + 

na^  e 

de 

dr 

de 

V  dr  J 
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dAC\ 
d?  ) 


de 


^  dAC\ 

de  dt  \  ^  ) 


di 

dt 


dVt 

dt 


COS  i 

'd{-An) 

na^  (1  —  e^)^A  sin  i 

du 

d AC\  dg 

$  + 

• 

V 

dr  J  du 

1 

'd{-An) 

na^  (1  —  e'^yA  sin  i 

dQ 

$  +  — ^ 

\  dr 

1 

'di-AH) 

na^  (1  —  sin  i 

di 

dAC  ^  dAC\ 
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df  d  U  dAC\ 
doj  dt  \  d?  ) 

dAC  dAC\ 

dr  dil  \  dr  ) 

dg  dfdf^  d  AC  \ 
(9f2  (9f2  (it  \  ,9f  y 


dAC  d  ,9A£\ 
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.9A£\ 
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^  dAC\ 
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1  -  \d{-An) 

no?  e  de 


dAC  ^  ,9A£\ 

dr  (9e  \  dr  ) 


dAC\ 

dr  / 
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^  dAC\ 

de  dt  \  (9F  / 


^  r  (9  (  -  AH) 
na  da 


dAC  d 
dr  da  y 


.9A£\ 
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,9A£\ 

~W  J 
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^  !L($  '9Ai:\ 
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(68) 


(69) 


(70) 
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Similarly,  the  gauge-invariant  Delaunay-type  system  can  be  written  down  as: 


dL  di-^n)  dAL  d  ( ^  dAC\  d AC\  dg  d?  d  d AC\ 

It  ~  ^  ^  ^ ~  ^  ^  j  ^  ^  ^  ^  ~W  j  ’ 


^  ^  dj-AH)  dAC  d  ^  dAC\ 

dt  dL  Qy  dL  y  Qy  j  \  dv  j  dL  dL  dt  y  Qy  j 

^  _  d  (  -  AH)  _  dAC  -  dAC\  _  /  -  dAC\  ^  ^  ±(^  dAC\ 

dt  doj  dr  doj  \  dC  J  y  dC  /  du  doj  dt  \  dC  ) 

duj  di  —  ALL)  dAC  d  d AC\  d AC\  dg  dfdf^  d AC\ 

*  "  +  +  i^  +  ^jaG  +  acSr  +  ^j 

^  _  d  ( -  AH)  _  dAC  -  dAC\  _  /  -  dAC\  ^  dAC\ 

dt  dVl  dtv  aS7  y  dtv  j  \  dr  j  ai2  aS7  dt  y  Qf  j 

dn  _  _a(-AW)  BAC  8  /-  aAC\  /-  dAC\  dg  dt  d  8AC\ 
dt  dH  dr  dH  y  dr  /  y  dr  /  OH  dH  dt  y  dr  / 

where 

L  =  /rd^ad^,  G  =  (l  -  e^)''\  H  =  (l  -  e^f'Los  i  (77) 


(72) 

,  (73) 
.  (74) 

(75) 
^  (76) 


— * 

and  the  symbols  $,  f,  g  denote  the  functional  dependencies  of  the  gauge,  position  and  velocity 
upon  the  Delaunay,  not  Keplerian  elements,  and  therefore  these  are  functions  different  from 
f,  g  used  in  dnn-izni)  where  they  stood  for  the  dependencies  upon  the  Kepler  elements. 
(In  Efroimsky  (2002)  the  dependencies  $,  f,  g  upon  the  Delaunay  variables  were  equipped 
with  tilde,  to  distinguish  them  from  the  dependencies  upon  the  Kepler  coordinates.) 

The  above  equations  do  not  merely  repeat  those  derived  earlier  in  Efroimsky  (2002,  2003), 
but  generalise  them  to  the  case  of  a  perturbation  AC  which  is  both  position-  and  velocity- 
dependent.  For  this  reason,  our  gauge-invariant  equations  can  be  employed  not  only  in  an 
inertial  frame  but  also  in  a  wobbling  one. 

To  employ  the  gauge-invariant  equations  in  analytical  calculations  is  a  delicate  task:  one 
should  always  keep  in  mind  that,  in  case  $  is  chosen  to  depend  not  only  upon  time  but  also 
upon  the  ’’constants”  (but  not  upon  their  derivatives),  the  right-hand  sides  of  these  equation 
will  implicitly  contain  the  hrst  derivatives  dCi/dt  ,  and  one  will  have  to  move  them  to  the 
left-hand  sides  (much  like  in  the  transition  from  (ll()|l  to  m)- 


Appendix  2: 

The  Hamilton- Jacobi  method  in  celestial  mechanics 


The  Jacobi  equation  (P|)  is  a  PDE  of  the  hrst  order,  in  (N-l-1)  variables  (g^,  t) ,  and  its 
complete  integral  lK(g,  Q,  t)  will  depend  upon  iV  -f  1  constants  a„  (Jeffreys  and  Jeffreys 
1972,  Courant  and  Hilbert  1989).  One  of  these  constants,  a^+i ,  will  be  additive,  because 
W  enters  the  above  equation  only  through  its  derivatives.  Since  both  Hamiltonians  are,  too. 
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defined  up  to  some  constant  /  ,  then  the  solution  to  ()d4|l  must  contain  that  constant  multiplied 
by  the  time: 


W(^q,  Oi,  ... ,  t)  —  hh(5',  cii,  ... ,  cini  t)  (t  to)  /(ui,  ... ,  cin) 


—  W [q,  Oi,  ... ,  (Iff,  t)  t  /(oi,  ... ,  Ojv)  fljv+i 


(78) 


where  the  hducial  epoch  is  connected  to  the  constants  through  to  =  —  ajv+i//  ,  and  the 
function  W  depends  upon  N  constants  only.  As  the  total  number  of  independent  adjustable 
parameters  is  N  +  1 ,  the  constant  /  cannot  be  independent  but  must  rather  be  a  function 
of  oi,  ... ,  a^,  ajv+i  .  Since  we  agreed  that  the  constant  a^+i  is  additive  and  shows  itself 
nowhere  else,  it  will  be  sufficient  to  consider  /  as  a  function  of  the  rest  N  parameters  only. 
(In  principle,  it  is  technically  possible  to  involve  the  constant  Ojv+i  ,  i.e.,  the  reference  epoch, 
into  the  mutual  transformations  between  the  other  constants.  However,  in  the  applications 
that  we  shall  consider,  we  shall  encounter  only  equations  autonomous  in  time,  and  so  there 
will  be  no  need  to  treat  Ojv+i  as  a  parameter  to  vary.  Hence,  in  what  follows  we  shall  simply 
ignore  its  existence.)  The  new  function  W  obeys  the  simplihed  Jacobi  equation 


/  dW{q,  ai,  a^,t) 

"A - Ff - 


+ 


dW{q,  ai,  ... ,  Ojv,  t) 
dt 


/(ai,  ...,  a^)  +  n*  (79) 


As  agreed  above,  is  a  constant.  Hence,  we  can  state  about  this  constant  all  the  same  as 
about  the  constant  / :  since  the  integral  W  can  contain  no  more  than  A^  +  1  adjustable 
parameters  oi,  ... ,  Ojv,  Ojv+i  ,  and  since  we  ignore  the  existence  of  a^v+i ,  the  constant  H* 
must  be  a  function  of  the  remaining  N  parameters:  Ti*  =  7-f*(ai,  ... ,  a  if). 

Now,  in  case  7i  depends  only  upon  (g,  p)  and  lacks  an  explicit  time  dependence,  then 
so  will  W ;  and  the  above  equation  will  very  considerably  simplify: 


(  dW{q,  ai,  ,  aif)\  ^  ,  n.j*r  \  (on\ 

n  I  g  ,  — I  =  /(ai,  ... ,  Ojv)  +  H  (ai,  ... ,  a^)  ,  (80) 

where  we  deliberately  avoided  absorbing  the  constant  Hamiltonian  Ti*  into  the  function  / . 

Whenever  the  integral  W  can  be  found  explicitly,  the  constants  (oi,  ... ,  Ojv)  can  be 
identihed  with  the  new  coordinates  Q  ,  whereafter  the  new  momenta  will  be  calculated  through 
P  =  —  dW/dQ  .  In  the  special  case  of  zero  Ti*  ,  the  new  momenta  become  constants,  because 
they  obey  the  canonical  equations  with  a  vanishing  Hamiltonian.  In  the  case  where  7Y*  is  a 
nonzero  constant,  it  must,  as  explained  above,  be  a  function  of  all  or  some  of  the  independent 

parameters  (oi,  ... ,  a^v)  ,  and  therefore,  all  or  some  of  the  new  momenta  P  will  be  evolving 

in  time. 

Since  it  is  sufficient  to  hnd  only  one  solution  to  the  Jacobi  equation,  one  can  seek  it  by 
means  of  the  variable-separation  method:  equation  (IHIHl  will  solve  in  the  special  case  when  the 
generating  function  (|7H|)  is  separable: 


N  _ 

IH(gi,  ... ,  gjv,  ai,  ... ,  a^)  =  ^  W  (g^ ,  Oi,  ... ,  a^v)  (81) 

i=l 

This  theory  works  very  well  in  application  to  the  unperturbed  (two-body)  problem  of 
celestial  mechanics,  a  problem  that  is  simple  due  to  its  mathematical  equivalence  to  the  grav¬ 
itationally  bound  motion  of  a  reduced  mass  'rripianet'^sunl {^planet  +  ^sun)  about  a  hxed 
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centre  of  mass  rripianet  +  •  If  one  begins  with  the  (reduced)  two-body  Hamiltonian  in 

the  spherical  coordinates 


qi  =  r  ,  =  (j)  ^  =  0  (82) 

(where  x  =  r  cos(j)  cos  6  ,  y  =  r  cos0  sin6'  ,  z  =  r  sin0),  then  the  expression  for  Lagrangian, 


^  ^  -  n  =  ^  {qif  +  ^  {qif  {q2f  +  ]-  {qif  {q^f  cos^  ga  +  — 

/  /  /  qi 


will  yield  the  following  formulae  for  the  momenta: 


dC 

Pi  =  ^  =  qi 

oqi 


dC  2  . 

P2  =  ^  =  qiq2 

dq2 


dC  2  •  2 

P3  =  ^  =  qiqs  cos  ga 
dqs 


whence  the  initial  Hamiltonian  will  read: 


w  =  Ep9  -  £  =  Ip?  +  ^p1  + 

Then  the  Hamilton-Jacobi  equation  (30)  will  look  like  this: 


2  P 

Ps - 

qi 


(83) 


(84) 


(85) 
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'dW' 

dqi. 


1  fdW' 


+ 


1  fdWy  _  y  _  dW 

2ql  cos^  ga  I  dq^,  I  gi  dt 


-  n*  =  0 


(86) 


while  the  auxiliary  function  W  dehned  through  dZHl)  will  obey 


1 
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'dW' 

.9qi, 


+ 


1  fdW' 


2gf  \dq2 


+ 


1  fdW\  _ 

2ql  cos2  ga  I  dqs  I  gi 


-  f  -  n*  =  0 


(87) 


A  lengthy  but  elementary  calculation  (presented,  with  some  inessential  variations,  in  Plummer 
1918,  Smart  1953,  Pollard  1966,  Kovalevsky  1967,  Stiefel  and  Scheifele  1971,  and  many  other 
books)  shows  that,  for  a  constant  H*  and  in  the  ansatz  (IHT|.  the  integral  of  (IHn|l  takes  the 
form: 


W  —  Wi  (gi,  Oi,  Oa,  03)  -|-  W2  (ga,  Oi,  CI2,  ^3)  +  kPs  (o's,  02,  03)  — 

(88) 

ei  (2{f  +  n*)  +  —  -  dqi  +  f  €2  (aj - ^  dga  +  /  03  dq^ 

to)  \  qi  qf  J  Jo  \  cos^  ga  J  Jo 

where  the  epoch  and  factors  ei,  ea  may  be  taken  as  in  Kovalevsky  (1967):  time  to  is  the  instant 
of  perigee  passage;  factor  ei  is  chosen  to  be  -|-  1  when  gi  =  r  is  increasing,  and  is  —  1  when  r 
is  decreasing;  factor  ea  is  -|-  1  when  ga  =  0  is  increasing,  and  is  —  1  otherwise.  This  way  the 
quantities  under  the  first  and  second  integration  signs  have  continuous  derivatives.  To  draw 
conclusions,  in  the  two-body  case  we  have  a  transformation-generating  function 


H^  =  IK  +  t /(ai,  ...,  a^)  = 


2(/  -f  n*)  + 


2/i 

qi 


-i]  dq,  + 

qi) 


(4  \ 

cos2  ga  J 


1/2 

dga  + 


(89) 


as  dqs  +  tf 


21 


whose  time-independent  component  W  enters  equation  iOD.  The  hrst  integration  in  (18911 
contains  the  functions  f(ai,  ,  a^)  and  ... ,  a^)  ,  so  that  in  the  end  of  the  day  W 

depends  on  the  N  constants  Oi,  ,  Ojv  (not  to  mention  the  neglected  to  ,  i.e.,  the  a^+i  ). 

Different  authors  deal  differently  with  the  sum  (/  -|-  Ti*)  emerging  in  (j89|l .  Smart  (1953) 
and  Kovalevsky  (1967)  prefer  to  put 

/  =  0  ,  n*  =  ai  ,  ai  =  -/x/(2a)  ,  (90) 

whereupon  the  new  momentum  Pi  =  —dW/dQi  =  —dW/dai  becomes  time-dependent 
(and  turns  out  to  equal  —t  +  to.)  An  alternative  choice,  which,  in  our  opinion,  better  reflects 
the  advantages  of  the  Hamilton- Jacobi  theory,  is  furnished  by  Plummer  (1918): 

f  =  ai  ,  H*  =  0  ,  oi  =  ^/J^d  .  (91) 

This  entails  the  following  correspondence  between  the  new  canonical  variables  (the  Delaunay 
elements)  and  the  Keplerian  orbital  coordinates: 


Qi  =  tti  = 

y/JTd 

Pi  =  -  Mo 

Q2  =  0,2  = 

(1  -  e2) 

P2  =  —  UJ  , 

(92) 

Q3  =  0,3  = 

V/i  a  (1  —  e2)  cos  i  , 

P3  =  -  G 

Everywhere  in  this  paper  we  follow  the  convention  (EH)  and  denote  the  above  variables  Qi,  Q2,  Qs 
by  L,  G,  H,  correspondingly  (as  is  normally  done  in  the  astronomical  literature) 
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